Nonlinear Relaxation Dynamics in Elastic Networks and Design Principles of 

Molecular Machines 
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Analyzing nonlinear conformational relaxation dynamics in elastic networks corresponding to two 
classical motor proteins, we find that they respond by well-defined internal mechanical motions to 
various initial deformations and that these motions are robust against external perturbations. We 
show that this behavior is not characteristic for random elastic networks. However, special network 
architectures with such properties can be designed by evolutionary optimization methods. Using 
them, an example of an artificial elastic network, operating as a cyclic machine powered by ligand 
binding, is constructed. 



INTRODUCTION 

Understanding design principles of single-molecule ma- 
chines is a major challenge. Experimental and theoreti- 
cal studies of proteins, acting as motors [l], 0, H, 0, @, 
ion pumps B, H, B or channels @, Q , and enzymes 




show that their operation involves 
functional conformational motions (see [H}). Such mo- 
tions are slow and cannot therefore be reproduced by full 
molecular dynamics simulations. Within the last decade, 
approximate descriptions based on elastic network mod- 
els of proteins have been developed H, 17 , 18 , TQ, 2oL 21|. 
In this approach, structural elements of a protein are 
viewed as identical point particles, with two particles 
connected by an elastic string if the respective elements 
lie close enough in the native state of the considered 
protein. Thus, a network of elastic connections corre- 
sponding to a protein is constructed. So far, the atten- 
tion has been focused on linear dynamics of elastic net- 
works, characterized in terms of their normal vibrational 
modes. It has been found that ligand-induced confor- 
mational changes in many proteins agree with the pat- 
terns of atomic disp lacements in their slowest vibrational 
modes [H H, Q (see also [H, 53) > even though 



nonlinear elastic effects must become important for large 
deviations from the equilibrium [28|, [22]. The focus of 
this article is on nonlinear relaxation phenomena in elas- 
tic networks seen as complex dynamical systems. 

Generally, a machine is a mechanical device that per- 
forms ordered internal motions which are robust against 
external perturbations. In machines representing single 
molecules, energy is typically supplied in discrete por- 
tions, through individual reaction events. Therefore, 
their cycles consist of the processes of conformational re- 
laxation that follow after energetic excitations. For a ro- 
bust machine operation, special nonlinear relaxation dy- 
namics is required. We expect that, starting from a broad 
range of initial deformations, such dynamical systems 
would return to the same final equilibrium state. More- 
over, the relaxation would proceed along a well-defined 
trajectory (or a low-dimensional manifold), rapidly ap- 
proached starting from different initial states and robust 



against external perturbations. These attractive relax- 
ation trajectories would define internal mechanical mo- 
tions of the machine inside its operation cycle. 

This special conformational relaxation dynamics has 
been confirmed in our study of the elastic networks of 
two protein motors (Fi-ATPase and myosin). On the 
other hand, our control investigation of random elastic 
networks has shown that relaxation patterns in random 
elastic networks are typically complex and qualitatively 
different from those of protein motors. Actual proteins 
with specific architectures allowing robust machine op- 
eration may have developed through a natural biological 
evolution, with the selection favoring such special dy- 
namical properties. In a model study, we have demon- 
strated that artificial elastic network architectures pos- 
sessing machine-like properties can be designed by run- 
ning an evolutionary computer optimization process. Fi- 
nally, an example of an artificially designed elastic net- 
work that operates like a machine powered by ligand 
binding has been constructed. 



Elastic network models 

The considered elastic networks consist of a set of N 
identical material particles (nodes) connected by identi- 
cal elastic strings (links). A network is specified by indi- 
cating equilibrium positions of all particles. Two parti- 
cles are connected by a string if the equilibrium distance 
between them is sufficiently small. The elastic forces, 
acting on the particles, obey Hooke's law and depend 
only on the change in the distances between them. In 
the overdamped limit [lq . the velocity of a particle is 
proportional to the sum of elastic forces applied to it. If 
R^ are equilibrium positions of the particles and Ri (t) 
are their actual coordinates, the dynamics is described 
by equations 
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FIG. 1: Elastic network of the single /3-subunit of the molecular motor Fi-ATPase (a) and the set of relaxation trajectories 

for this network (b). Links are colored according to their relative deformations Py'l m the motion corresponding to the 

slowest normal mode. Each of 100 trajectories starts from a different initial conformation (see Methods). In this case, all the 
trajectories converge to the original equilibrium state (blue dot). The labels (1,2,3) are attached at Ile390, Argl91 and Gly54, 
respectively. 



where A is the adjacency matrix, with the elements Ay = 



1, if 
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< Iq, and Ai 



otherwise. The 



dependence on the stiffness constant of the strings and 
the viscous friction coefficient of the particles is removed 
here by an appropriate rescaling of time. 

The dynamics of elastic networks is nonlinear, because 
distances |R, — Rj| are nonlinear functions of the coor- 
dinates Ri and Rj. Close to the equilibrium, equations 
of motion can however be linearized, yielding 
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for small deviations r, = R, R 
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These equations 



can be written as f j = — ■ ^-ij v ji where A is a 3 AT x 
3iV linearization matrix. In the linear approximation, 
relaxation motion is described by a sum of independent 
exponentially decaying normal modes 
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with A a and representing nonzero eigenvalues and 
the respective eigenvectors of the matrix A. It should 
be noted that the same eigenvalues determine vibration 
frequencies of the network, w a ~ \fK^, and the vibra- 
tional normal modes are the same. Long-time relaxation 
is dominated by soft modes with small eigenvalues. 



RESULTS 
Two motor proteins 

As an example, Fig. 1 shows the elastic network of 
the single /3-subunit of the molecular motor Fi-ATPase 



(Protein Data Bank ID: 1H8H, chain E) [30|. Each node 
corresponds to a residue in this protein (the total num- 
ber of nodes is N — 466). In Fig. lb, a set of confor- 
mational relaxation trajectories, obtained by numerical 
integration of the nonlinear elastic model (see Methods, 
Eq. (1)) of this macromolecule, is displayed. To track 
conformational motions, three network nodes (1, 2 and 
3) were chosen and pair distances U12, U13 and U23 were 
determined during the relaxation process. Thus, each 
conformational motion was represented by a trajectory in 
a three-dimensional space, where coordinates were nor- 
malized deviations Au^/u^ ' 



from the equilibrium pair 



distances 



,(°) 



Auij = Uij — ufP)- Each trajectory begins from a dif- 
ferent initial conformation obtained by applying random 
static forces to all network nodes (sec Methods). Trajec- 
tories starting from various initial conditions soon con- 
verge to a well-defined relaxation path leading to the 
equilibrium state. This path corresponds to a slow mo- 
tion of the network group, including label 1, with respect 
to the rest of the molecule. The protein is soft along such 
a path: by applying static forces of the same magnitude, 
one can stretch it by 30% along the path, as compared to 
the length changes of only a few percent when the forces 
were applied in other directions. 

Relaxation behavior in the nonlinear elastic network 
(N = 793) of another classical molecular motor, myosin 
(single heavy chain; PDB ID: 1KK8, chain A) [31[, is 
displayed in Fig. 2. In contrast to the /3-subunit of 
Fi-ATPase, this elastic network possesses an attractive 
two-dimensional manifold (a plane). The network is ex- 
tremely stiff for deformations in the directions orthogo- 
nal to this plane. By applying static forces of the same 
magnitude, one can only induce relative deformations of 
about 10 -3 along such orthogonal directions, as com- 
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FIG. 2: Elastic network of the single heavy chain of myosin (a) and the set of relaxation trajectories for this network (b,c,d). 
Panels (b) and (c) are viewed from different angles; each of 100 trajectories starts from a different initial conformation (see 
Methods). All the trajectories converge to the original equilibrium state (blue dot). Panel (d) shows a trajectory (shown 
by green curve in (b) and (c)) with labels of passage time. The labels (1,2,3) are attached at Leu836, Asp63 and Glu370, 
respectively. 
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FIG. 3: Eigenvalue spectra for the networks of Fi-ATPase 
and myosin, normalized to the lowest nonzero eigenvalue Ai. 
Ai = 1.09 x 10 -2 for Fi-ATPase and 2.81 x 10~ 5 for myosin. 



pared to the relative deformations of the order of 10 _1 
for the directions within the plane. To characterize the 
temporal course of relaxation, one relaxation trajectory 
(displayed by green color in Fig. 2 (b,c)) and subsequent 
positions at different time moments along the trajectory 
are indicated in Fig. 2d. The trajectory rapidly reaches 
the plane and then the relaxation motion becomes much 
slower (with the characteristic time scales larger by a 
factor 10 to 100). A similar behavior is characteristic for 
other relaxation trajectories, starting from different ini- 
tial conditions. All recorded trajectories returned to the 
equilibrium state and no metastable states were encoun- 
tered starting from the chosen initial conditions. 

Nonlinear effects were essential in the relaxation dy- 



namics starting from large arbitrary initial deformations 
considered here. Remarkably, the observed relaxation 
patterns are nonetheless qualitatively in agreement with 
the normal mode analysis. Both motors possess a group 
of soft modes separated by a gap from the rest of the spec- 
trum (eigenvalue spectra of elastic networks of these pro- 
teins are shown in Fig. 3). The two soft modes of myosin 
define the attractive plane seen in the relaxation pattern 
of its elastic network (Fig. 2). The elastic network of 
Fi-ATPase has three soft modes. However, one of the 
soft modes has the relaxation rate which is smaller than 
the other two modes. Therefore, the pattern of relax- 
ation trajectories looks here like a thick one-dimensional 
bundle. Specific ligand-induced conformational changes 
in Fi-ATPase and myosin were previously shown to have 
strong overlaps with patterns of deformation in slow vi- 
brational modes 1241 1. 



Random and designed elastic networks 

A control study of nonlinear relaxation phenomena 
in random elastic networks was performed. Such net- 
works were obtained by taking a relatively short chain of 
N = 64 and randomly folding it in absence of energetic 
interactions (see Methods). After that, all particles sepa- 
rated by short enough distances were connected by iden- 
tical elastic links. Figure 4 shows relaxation patterns for 
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FIG. 4: Relaxation trajectories for two random elastic net- 
works, in the plane (Aun/ii$ , Auiz/ii]^) of normalized dis- 
tance deviations. Each of 100 trajectories starts from a differ- 
ent initial conformation. Blue dots indicate stationary states 
reached; the original equilibrium state is Aui2 = AU13 = 0. 
These networks have no internal rotational modes. 



two typical random elastic networks (using the same pro- 
cedure for generation of initial conditions and for tracking 
of conformational relaxation as in Figs. 1 and 2). Relax- 
ation patterns in random networks are clearly different 
from those of the considered motor proteins. Random 
networks possess many (meta)stable stationary states, 
with relaxation trajectories often ending in one of them 
instead of going back to the equilibrium conformation. 
The linear normal mode description holds in such net- 
works only in close proximity of the equilibrium state. 

Thus, elastic networks of motor proteins are spe- 
cial. Their equilibrium conformation has a big attraction 
basin. Starting from an arbitrary initial deformation, re- 
laxation dynamics is soon reduced to a low-dimensional 
attractive manifold. Within its large part, the dynam- 
ics is approximately linear and determined by a few soft 
modes. Proteins with such special dynamical properties, 
essential for their functions, may have emerged through a 
biological evolution. Below, we show that artificial elas- 
tic networks with similar properties can be constructed 
by running an evolutionary optimization process based 
on a variant of the Metropolis algorithm. 

The evolutionary optimization technique is described 
in Methods. For each network, spectral gap g = 
log 10 (A2/A1) is defined as the logarithm of the ratio be- 
tween the relaxation rates A2 and Ai of its two slowest 
normal modes. If a substantial gap is present, the slowest 



mode has the relaxation rate which is much smaller than 
that of the other modes; therefore, the long-time relax- 
ation in the linear regime would be dominated by this 
soft normal mode. The employed evolutionary optimiza- 
tion process maximizes the spectral gap g of the evolving 
networks. Beginning with an initial random network, we 
applied structural mutations and determined the differ- 
ence Ag = g' — g of the gaps before and after a mutation. 
If the gap was increased (Ag > 0), the mutation was al- 
ways accepted. If Ag < 0, the mutation was accepted 
with the probability P = exp (Ag/9) where 9 is the ef- 
fective optimization temperature. This procedure was 
applied iteratively, generating an evolution started from 



an initial network 35 1 



Networks with soft normal modes and large spectral 
gaps were thus constructed. A typical network with a 
large gap and its relaxation pattern are shown in Fig. 
5 (a,b) (for other examples, see Supplementary Fig. 1). 
The presence of a large gap has a strong effect on the 
nonlinear relaxation properties in such systems. They 
possess well-defined long paths with slow conformational 
motion leading to the equilibrium state. There are only 
a few (meta) stable states and these states usually lie on 
an attractive relaxation path, so that a small barrier is 
encountered when moving along it. The opposite behav- 
ior with complex relaxation patterns and a high num- 
ber of (meta) stable conformations was found for a set of 
"failed" networks where gaps were small and could not be 
significantly increased through the evolution (see Supple- 
mentary Fig. 2). Our analysis shows that the designed 
networks can be viewed as consisting of rigid blocks con- 
nected by soft joints; they are able to perform some large 
conformational changes accompanied by only small local 
deformations (see Supplementary Fig. 3). 

Spectral gaps of the networks, designed by using this 
optimization method, are much larger than those char- 
acteristic for motor proteins (cf. Fig. 3). To improve 
the agreement, additional "reverse" evolution was sub- 
sequently applied to the designed networks, with the se- 
lection pressure aimed to decrease the gap (see Meth- 
ods). While the gap was rapidly reduced, the global 
relaxation pattern was changing only much more slowly 
with structural mutations and retained characteristic fea- 
tures of the networks with large gaps. Figure 5c displays 
the network, obtained by applying such reverse evolution 
(with only 5 subsequent mutations) to the original net- 
work shown in Fig. 5a. Although the spectral gap has 
been reduced from 9.53 to 1.22, the principal structure 
of the network is retained, with the mutations mostly 
affecting only the hinge region. The relaxation pattern 
of the constructed network (Fig. 5d) reveals an attrac- 
tive path leading to the equilibrium state (with another 
stable state lying on it). Remarkably, the linear normal 
mode description of relaxation dynamics holds in such 
networks within a much larger domain around the equi- 
librium state. We conclude that, by running an evolu- 
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FIG. 5: Examples of designed elastic networks, (a) A typical network with a large gap (g — 9.53) and (b) its set of 100 
relaxation trajectories; (c) a network with a smaller gap (g — 1.22) obtained by reverse evolution from the network in panel 

(a), and (d) its set of 100 relaxation trajectories. Network links are colored according to their relative deformations \p2^\ m 
the motion corresponding to the slowest normal mode; thick lines indicate the backbone chain. Each trajectory starts from a 
different initial conformation; blue dots indicate stationary states reached. 



tionary optimization process, artificial elastic networks 
approaching conformational relaxation properties of real 
motor proteins can be constructed. 



An artificial machine network 

Using such designed networks, we proceed to construct 
nonlinear elastic systems which can be viewed as pro- 
totypes of a machine powered by ligand binding. The 
network shown in Fig. 6 performs cyclic hinge motions, 
caused by binding and detachment of a ligand. To ob- 
tain it, an initial network with two distinct parts (dense 
clusters), which were only loosely connected, was pre- 
pared. This initial network was characterized by a small 
spectral gap. By running evolutionary optimization, the 
architecture of the network was modified, so that it de- 
veloped a large spectral gap while maintaining its special 
structure. The equilibrium conformation of the finally 
obtained network is shown in the lower right panel in 
Fig. 6. The cycle began with binding of a ligand (snap- 
shot A, t = 0). The ligand was an additional particle that 
could form elastic connections with the existing network 
nodes. Binding of a ligand was modeled by placing the 
particle at a location chosen in the hinge region and cre- 
ating elastic links with the three near nodes (for details, 
see Methods). When such new links were introduced, 



they were in a deformed (i.e., stretched) state and a cer- 
tain amount of energy was thus supplied to the system. 
Therefore, the network-ligand complex was not initially 
in the equilibrium state and conformational relaxation to 
the equilibrium had to occur. Within a short time, de- 
formations spread over a group of near elastic links (B, 
t = 200), producing a small change in the network config- 
uration. After that, a slow hinge motion of the network- 
ligand complex took place (C, t = 4000) until the equi- 
librium state of the complex was reached (D, t = 20000). 
fn this state, the ligand was removed. At this moment, 
the elastic network was in a configuration far from equi- 
librium and its relaxation set on. Within a short time, 
the network reached the path (E, t = 20200) along which 
its slow ordered relaxation proceeded (F, t = 30000), un- 
til the equilibrium conformation (A) was reached. Then, 
a new ligand could bind at the network, and the cycle 
was repeated. For visualization of conformational mo- 
tions inside the network cycle, see Supplementary Video 
1. 

To monitor conformational motions, three labels were 
attached to the network and distances U12 and U13 be- 
tween them were recorded. The lower left panel in Fig. 6 
shows trajectories in the plane (Awi2,Aui3) which corre- 
spond to the operation of this prototype machine in the 
presence of noise and with stochastic binding of ligands. 
As seen in Fig. 6, the trajectories of the forward and back 
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FIG. 6: The prototype of a molecular machine: an elastic network performing regular cyclic hinge motions caused by repeated 
binding and detachment of ligand particles. A series of snapshots (A-F), corresponding to different cycle phases, is displayed. 
In the snapshots, the stretched links are shown in light blue (0.01 < Auy < 0.1) and in blue (0.1 < Auij), the compressed 
links in orange (—0.1 < Awy < —0.01) and in red (Auy < —0.1). The ligand is shown as a red dot. The lower left panel shows 
trajectories of the network for 10 subsequent cycles in the plane (Aui2,Aui3) of distance deviations between the three labeled 
nodes, indicated in the right panel. Dots A to F along the trajectory correspond to the network conformations shown in the 
snapshots A to F above. 



hinge motions are different, but both of them are well 
defined. The applied noise (see Methods) is only weakly 
perturbing them. This means that both the free network 
and the network-ligand complex have narrow valleys with 
steep walls, leading to their respective equilibrium states. 
Binding of a ligand is possible in a relatively broad inter- 
val of conformations near the equilibrium and therefore 
there is a dispersion in the transitions from the upper to 
the lower branch of the cycle. The operation of this ma- 
chine is essentially nonlinear and unharmonic effects are 
strong. Characteristic rates of slow relaxation motions in 
this network differ by a factor between 10 and 1000 from 
the rate of the fast conformational transition following 
the ligand binding. 



DISCUSSION 

In conclusion, we have shown that motor proteins pos- 
sess unique dynamical properties, intrinsically related to 
their functioning as machines. We have also demon- 
strated that artificial elastic networks with similar prop- 
erties can be constructed by evolutionary optimization 
methods. To verify these theoretical predictions, special 
single-molecule experiments monitoring conformational 
relaxation after arbitrary initial deformations can be per- 
formed. An example of an elastic machine-like network 
powered by ligand binding is presented. Using such de- 
signed networks, fundamentals of molecular machine op- 
eration, including energetic aspects, the role of thermal 
noise and hydrodynamic interactions, can be discussed 
[3^ |. Comparing the behavior of such designed networks 
with that of real molecular machines, better understand- 
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ing of what is general and what is specific for a particular 
protein can be gained. Moreover, our analysis provides a 
systematic approach for the design of proteins with pre- 
scribed (programmed) robust conformational motions. 
Not only proteins, but also atomic clusters can be de- 
scribed by elastic network models (3^|- Similar meth- 
ods can further be used for engineering of non-protein 
machine-like nanodevices (see [34|). 



METHODS 

The elastic networks for protein structures in Figs. 1 
and 2 are constructed from the structural data of Fx- 
ATPase in the ATP-analog state (Protein Data Bank ID: 
1H8H) chain E [3(| and myosin in the ATP-analog state 
(PDB ID: 1KK8) chain A |3l| with a cutoff distance l = 
10 A. We place a node at the position of each a-carbon 
atom, and connect nodes lying within the cutoff distance 
by a link. Stiffness constants of all links and friction 
coefficients for all particles are equal. 

To generate random networks, a chain of N nodes 
is taken and folded randomly in the three-dimensional 
space. Each next node is positioned at random, in such 
a way that its distance I from the preceding node satis- 



node should be separated by at least the distance Z m j n 
from all previous nodes. Examining the folded chain, all 
pairs of nodes with distances less than lg are noticed 
('o > 'max) and connected by elastic strings (hence, not 
only all neighbors in the string get connected, but also 
those nodes which come by chance close one to another 
in the folded conformation). 

The iterative optimization process consists of a se- 
quence of structural mutations followed by selection. To 
perform a mutation, a node is taken at random and its 
equilibrium position is changed. The new equilibrium 
position is chosen with equal probability within a sphere 
of radius d around the old equilibrium position. To pre- 
serve the backbone chain, we additionally require that, 
after a mutation, the distances of the node from its both 
neighbours in the chain lie within an interval from l m [ n to 
'max! other pair distances should not be shorter than Z m ; n 
either. The new graph of connections after a mutation is 
constructed by examining the pair distances between all 
nodes and linking the nodes separated at equilibrium by 
a distance shorter than Iq. 

Sometimes, networks allowing internal free rotations 
(and, thus, additional zero eigenvalues) may be obtained. 
To distinguish zero eigenvalues numerically, we adopt a 
numerical cutoff § A = 10~ 12 and assume the eigenvalues 
less than SX to be zero. For the networks without inter- 
nal rotation modes, that is, if the number N nz of nonzero 
modes is equal to 37V — 6, we proceed in each iteration 
step as described in the main text. When such modes are 
present, a mutation is always accepted when it decreases 



the number of zero eigenvalues {N' zm < 
cepted with the probability P = exp [— (N' zm — N zm ) /$] 
otherwise. In this study, we consider networks with 
N = 64 nodes. The parameter values d — l m i n = 3.4, 
'max = 4.2, Zq = 8.0 and 9 = 0.1 are used. 

Most of the networks after optimization (97.3%) had 
a large gap g > 3, while such gaps were rare in the ini- 
tial networks (1.9%). Here, only networks without inter- 
nal rotations were counted (1511 random networks out 
of 30000 trials and 2346 selected networks out of 2500 
trials). 

In the "reverse" evolution used to generate special net- 
works with smaller gaps (such as shown in Fig. 5c), the 
same evolution algorithm is employed with the replace- 
ment of g by — g; the effective temperature is = 0.01 in 
these simulations. 

Around the equilibrium conformation, relative changes 
in the distances between nodes i and j in a relax- 



ation mode a are calculated as p\* — duij/dX a = e\" 

(0) / (0) , (a) (a) (a) 

u)//u)/, where e - " 



and u 



(o) 



(0) 



(a) 
ij 



y » 3 

Here, the eigenvec- 



tors are normalized as 



3 ' 

(see also Eq. ©). 

2 

— 1. In Supplementary 

Fig. 3, statistical distributions of relative changes p^ 
in the slowest relaxation mode (a = 1) in ensembles of 
elastic networks are shown. 

For trajectory visualizations, three labels (1,2,3) are 
attached to a network and distances Ui2, it 13 and U23 
between them are recorded. The labels are chosen in 



such a way that 



, the relative distance change in the 

slowest relaxation mode, is maximal between the nodes 
labeled as 1 and 2: then the third node is chosen for 



which 



(2) 
Pl3 



in the second slowest relaxation mode, is 

maximal between the nodes labeled as 1 and 3 (there are 
two choices and we have selected one of them) . 

To prepare initial conditions when relaxation patterns 
are considered, we apply randomly chosen static forces 
Fj lS taii C and wait for a certain time r s . These random 
forces are chosen in such a way that their total mag- 



nitude F s 



1/2 



is fixed. Then, we re- 



static] 

move the static forces and record the trajectory. For each 
network, 100 trajectories from different initial conditions 
are shown. In Fig. 1, F s — 10, t s = 3 x 10 4 ; in Fig. 2, 
F s = 0.1, t s = 10 5 ; and in Figs. 4 and 5, Supplementary 
Figs. 1 and 2, F s = 0.1 and t s = 10 6 . 

In Fig. 6 and in Supplementary Video 1, we initially 
place a ligand in the center of mass of the three ligand- 
binding nodes and introduce additional links to these 
nodes with the natural length of 1.7. At the initial mo- 
ment, the links are longer than their natural lengths (i.e., 
they are stretched) and thus attractive forces between the 
ligand and the binding nodes are acting. To include fluc- 
tuations, random time-dependent forces fj(t) of intensity 
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a have been added to the equations of motion of all par- 
ticles, i.e., Rj = F he ia S tic + fi(t) with (ii(t)) = and 
(f;(t) • fj(i/)) = 6aSijS(t - t'). We have a = 10~ 3 (upper 
panel in Fig. 6 and Supplementary Video 1) and 3 x 10 -3 
(lower left panel in Fig. 6). At t = 20000 (which is long 
enough for the network to relax to the steady state with 
a ligand, snapshot D), we remove the ligand and the ad- 
ditional links. In the ligand-free state, we have assumed 
that the network binds with a ligand again at a constant 
rate (probability per unit time) v even in noncquilibrium 
network configurations, provided that all distances Uij 
between the three ligand-binding nodes satisfy the condi- 
tions lAuiil = 



,(0) 



< e;, i.e., the conformation of 

the ligand-binding site is close enough to the initial one. 
Then, a ligand is placed again in the center of mass of 
the ligand-binding nodes and the next cycle starts. The 
parameter values v = 10~ 6 and e; = 0.01 are used. Co- 
ordinates of the nodes in the constructed example of the 
machine-like network and positions of the three binding 
nodes are given in Supporting Information Data 1. 
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Supplementary Figure 1 : Relaxation trajectories for networks with the largest gaps g (the 
best six cases) in 939 optimized networks. 100 trajectories for each network are shown 
in the plane (Au 12 , Au 13 ) of distance deviations between three selected nodes, 
normalized to the initial distances u 12 <°> and u 13 (°). Each trajectory starts from a different 
initial conformation obtained by applying random static forces to the network (in the 
same way as in Fig. 5; see Methods). Blue dots indicate final states; the original 
equilibrium state is Au 12 = Au 13 = 0. The gaps g are (a) 9.75, (b) 9.74, (c) 9.69, (d) 9.67, 
(e) 9.62, (f) 9.53. 




Supplementary Figure 2: Relaxation trajectories for networks with the smallest gaps g 
(the worst six cases) in 939 optimized networks. 1 00 trajectories for each network are 
shown in the same way as in Supplementary Fig. 1. Each trajectory starts from a 
different initial conformation obtained by applying random static forces to the network. 
Blue dots indicate final states; the original equilibrium state is Au 12 = Au 13 = 0. The gaps 
g are (a) 0.31, (b) 0.68, (c) 1.09, (d) 1.39, (e) 1.55, (f) 1.87. The relaxation processes 
are complex in these systems, and qualitatively different from those with a large gap 
shown in Supplementary Fig. 1 . 




Supplementary Figure 3: Distributions of relative changes |PijC)| in pair distances 
between the linked (red) and nonlinked (blue) nodes in the slowest normal mode (a) for 
the ensemble of 939 optimized networks and (b) the control ensemble of 939 networks. 
In the networks optimized for a large gap, a large fraction of all elastic links remains 
almost not deformed by the slowest normal mode motion. This indicates that such 
networks contain rigid blocks. Despite only very weak deformations inside the blocks, 
the constructed networks can show large shape changes. This is seen in the distribution 
of distance changes for the pairs of nodes not directly connected by elastic strings. 
Hence, the network functions as if it were constructed of rigid mechanical parts with 
flexible connections. The control ensemble is obtained by running evolution with 
mutations, but in absence of any selection on the gap g. 



